Quantum dynamics by the constrained adiabatic trajectory method 
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We develop the constrained adiabatic trajectory method (CATM) which allows one to solve the 
time-dependent Schrodinger equation constraining the dynamics to a single Floquet eigenstate, as 
if it were adiabatic. This constrained Floquet state (CFS) is determined from the Hamiltonian 
modified by an artificial time-dependent absorbing potential whose forms are derived according 
to the initial conditions. The main advantage of this technique for practical implementation is 
that the CFS is easy to determine even for large systems since its corresponding eigenvalue is well 
isolated from the others through its imaginary part. The properties and limitations of the CATM 
are explored through simple examples. 



I. INTRODUCTION 

Modern developments and applications of quantum 
mechanics often involve complex chemical and even bio- 
logical systems driven by laser fields (see for instance 
Solving (numerically) the time dependent Schrodinger 
equation (TDSE) for such time dependent systems be- 
comes then very time consuming and sometimes even 
impossible. Finding numerical simplifications is an ac- 
tive research. One can for instance mention the multi- 
configuration time-dependent Hartree (MCTDH) algo- 
rithm [2[. 

Techniques that lead to an efficient propagation of a 
time-dependent problem often involve the Floquet the- 
ory which allows one to incorporate fast oscillations of 
the external field (for instance such as the optical oscilla- 
tions of a laser field) in an enlarged Hilbert space [3| . For 
instance, it permits an adiabatic separation between the 
fast field oscillation dynamics and the slow time modula- 
tion of the field envelope (adiabatic Floquet theory [4|]). 
This Floquet technique can be used alternatively to treat 
the full time-dependence of the field, which is referred to 
as the (t,t r ) approach [f|. 

Relevant processes are most often expected to be de- 
scribed in a small subspace, often named active space, 
through effective Hamiltonians. One can mention in par- 
ticular the time-dependent wave operator theory (TD- 
WOT) as a tool to extract dynamical active spaces @. 

A few years ago Jolicard et al. Q have proposed the 
"Constrained Adiabatic Trajectory Method" (CATM) 
for solving the TDSE for a time-dependent potential. 
Since we use a quantum mechanical approach the tra- 
jectory studied in the CATM is not a classical one but 
rather a constrained path followed by the wavefunction as 
it develops in time in a composite Hilbert space which we 
describe below. We here investigate that method extend- 
ing it for an initial condition as a general superposition of 
states for a small system, and emphasizing its principal 
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novel feature, the use of a complex absorbing potential 
which is itself time-dependent. The usual approach is 
to propagate the wavefunction in small time steps, with 
the Hamiltonian considered as constant over each step 
d i]. The radically different approach of the CATM is 
to limit the time development to only one term in a Flo- 
quet expansion of the wavefunction, achieving this by a 
careful choice of the varying absorbing potential. The 
problem of integrating the TDSE then becomes that of 
finding one eigenvector of the system's Floquet Hamil- 
tonian. The method has some affinities with the (t, t') 
approach Q but represents a modification and improve- 
ment of it. The method finds the wavefunction at regular 
grid points throughout the interaction period, the prin- 
cipal requirement being to work with a sufficient number 
of points to describe the time-varying Hamiltonian and 
to allow the stable use of Fast Fourier transforms. 

In brief, the technique requires the wavefunction corre- 
sponding to the dynamics to connect to a single Floquet 
state, referred to as a constrained Floquet state (CFS), 
through the use of an artificial absorbing potential (or 
optical potential). The second role of the absorbing po- 
tential is to dilate the Floquet spectrum isolating well the 
eigenvalue corresponding to the CFS from the other ones. 
In practice one thus needs to find this CFS to determine 
the dynamics. 

In Section II, on the basis of Ref. Q , we summarize the 
technique CATM with its corresponding Floquet repre- 
sentation, and recall the result when the initial condition 
is a single eigenstate of the free system. In Section III, 
we extend the technique to a more general initial state, 
as a superposition of eigenstates of the free system. This 
is analyzed for a two-state system. A forthcoming work 
will treat the case of systems of higher dimension. In Sec- 
tion IV, we give an analytic treatment of the effect of the 
absorbing potential on the Floquet spectrum for a two- 
state model. The numerical limitations of the method 
and its accuracy are analyzed in section V. We study ex- 
amples with two- and three-level models which illustrate 
the dual role played by the optical potential in Section 
VI. Section VII is devoted to the conclusion. 
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II. THE CONSTRAINED ADIABATIC 
TRAJECTORY METHOD 

A. The Floquet representation 

We assume a system of Hamiltonian H (q, t) (where the 
quantum coordinates have been denoted by q) defined in 
a basis This Hamiltonian can be usually decom- 

posed as H (g, t) = H (q)+W(q, t) featuring a free system 
Ho(q) subjected to an external time dependent field cor- 
responding to the interaction potential W(q,t). In that 
case, {\j)} correspond to the states of the free system. 
We assume that the interaction potential W(q, t) acts on 
a duration t £ [0, T] referred to as the physical duration 
in the following. We define an extra time interval [T, T'] 
after the physical interaction time during which (i) we 
add an artificial time-dependent absorbing (or optical) 
potential V(q,t) satisfying V(q,0 <t<T) = V(q,T') 
. and (ii) we extend continuously the interaction such 
that W(q,T') — W(q, 0). This construction features a 
periodic Hamiltonian H(q,T') = H(q,0). 

We can thus define the corresponding Floquet Hamilto- 
nian (or quasi-energy operator) on the extended Hilbert 
space (product of the original Hilbert space, i.e. asso- 
ciated to the free system, by the space of T'-periodic 
functions) (ij: 

H F (q, t) = H (q) + W(q, t) + V(q, t) - ih^. (1) 

We consider the entire duration of the interac- 
tion+absorbing potential as a fundamental period T" 
(luq = 2n/T'), contrary to the traditional Floquet scheme 
in which T" is associated with the period of an external 
field (such as the optical period of a laser field). The 
Floquet states can be indexed with a double labelling 
j, n linked to a finite basis representation of the decou- 
pled parts of ([1]), i.e. to the free-system eigenstates 
(j \j)) an d to the operator —ihd t (corresponding here 
to a Fourier basis, n -H- \n) = \e~ mu)at )). Thus a complete 
basis is formed with the eigenstates {\Xj }U (q, t))} of Hp : 

Hp\X jt7l (q,t)) =Hu Xi J\j,n(q,t)). (2) 

Using the periodicity properties of the Floquet theorem 
{\Xj, n (q,t)) = \\ j ,o(q,t))e* nu <> t , u Xin = ^ +n Wo ), we 
can rigorously expand the solution of the time dependent 
Schrodinger equation (TDSE) with an initial limitation 
to the first Brillouin zone Q, 

|*(?,t)) = $> ij0 (g, 0))e-^.o t |A j , o (g,t)>- (3) 

3 

(Here we consider for simplicity only a bound spectrum, 
that can feature however imaginary parts; the exten- 
sion to a system with bound and continuous spectrum is 
in principle direct assuming a discretization of the con- 
tinua). Usually, a great number of | Aj ; o) is necessary to 
reconstruct \*£(q, t)}. An interesting practical application 
of Eq. ^ is the development of a very reduced number of 



Floquet vectors, and in the best case of only one, which 
is the key idea of the CATM. 

The method developed in this paper deals with the case 
of a single constrained Floquet state (CFS) and labeled 
I in the expansion ([3]). In this case, the CFS has to 
match, when projected at t = 0, with the initial boundary 
conditions required for the wavefunction ^{t = 0): 

(A i , (g,0)|*(g,0))=^, (4) 

i.e. 

\*(q,t))=e- i ^ t \\(q,t)), (5) 

where we have omitted in the latter the index I for sim- 
plicity: w A = u)xi, , |A(g, t)) = \Xi t0 (q,t)). 

We will show below that in practice we do not get the 
exact equality ([5]) but a proportionality through a well 
defined complex phase. 



B. Initial condition as an eigenstate of the free 
system 

Jolicard et al. 0] provided the matching with the ini- 
tial condition for an initial state equal to a single state 
|i) of the {\j}} basis, i.e. \^(q, 0)) = The connection 
between the Floquet eigenstate and the required initial 
state is made thanks to the addition of the absorbing 
potential V on the extra interval [T, T'] . Below we sum- 
marize this scheme and extend it in the next section to 
any required initial condition for the particular case of a 
two-state system. 

In order to satisfy Eq.([S]) (with a proportionality in- 
stead of the equality), it is sufficient to have the connec- 
tion at t = 0: 

|A(?,0))oc|i>. (6) 

We remark that \X(q,t)) is a Floquet vector of the ex- 
tended Hilbert space, but fixing t to a specific value leads 
to a component of this vector of dimension of the origi- 
nal Hilbert space. Eq. © suggests the use of the following 
form for the absorbing potential: 

V(t)=J2-iVoAt)\m (7) 

with 14pt(£) zero over [0,T] and positive over [T,T'}. As 
shown in [7|, provided that 

\J V opt {t)dt^ \$S(u x )\(T'-T) (8) 

we can be sure that all channel except \i) are absorbed 
and that Eq.Q is satisfied to a good approximation (as 
will be tested in section lyi) . 
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III. EXTENSION OF THE CATM TO A 
GENERAL INITIAL CONDITION: THE 
TWO-STATE CASE 

If we wish to work with CATM in the case of an initial 
condition as a state superposition of the free system, i.e. 



l*(Q)>=$>|j), 



(9) 



then simple forms as (|7J| no longer work. (From now 
on, we do not write explicitly the dependence on the q 
coordinates.) 

We provide below the relevant absorbing potential that 
should be used for a two-state system of Hamiltonian 



H(t) 



Ai(t) Q(t) 

n*(t) A 2 (t) 



(10) 



We consider the most general case with diagonal terms 
that are time dependent (due to Stark shifts of the 
states for instance) and complex (i.e. including their life- 
time). We assume that the coupling is in general 
different from zero only over the physical time interval 
[0, T]. During the extra time interval, the diagonal terms 
have to be continuously varied such that they recover 
their initial value in order to guarantee the periodicity: 
A J -(T') = A i (0). 

We show below that it is possible to treat any ini- 
tial condition by adding the following absorbing potential 
over the supplementary interval [T, T'] : 



/ 

V(t) = cae n?w)«' 

\ cie isr a X (o«' 

with V opt (t)>Q Vte}T,T'[ 
V opt (t)=0 ViG[0,T]. 



x {-iV opt (t)) (11) 



The operator 



n(t) = 



i/t T A 2 (t')it' 
i/ t T ' Ai (*')<«' 



(12) 



involved in this definition (1111) is a non-orthogonal (i.e. 
non self-adjoint) projector, i.e. II 2 = II, whose kernel is 
the initial state up to a phase correction: 



/ Cie *j; T Alto*' \ 
n(t Cie L , , = o. 



(13) 



For this case it is indeed possible to obtain the analytical 
asymptotic form of the Floquet eigenvector over the ex- 
tra interval [T, T'] , where the Hamiltonian contains just 
the free system Hamiltonian and the absorbing poten- 
tial. With the above definition and writing Floquet com- 
ponents (j'|A(i)) = Xj(t), one must solve on [T, T'\ the 



following system : 
d\i{t) 



Of. 



i(w A - A 1 (t))Ai(t) 



(14a) 



d\ 2 (t) V opt (t) c 2 e'/t T WW 



dt 



■ Cie if t T W)dt> 

V opt (t) 



Ai(t) 



-i(w x -A 2 (t)) A 2 (t) (14b) 



The first component follows an exponential law: 
Ai(t) = Ai(T)e i ( a, ^ t - T )-i'T A i (*')*'). This function 
can be introduced in the second equation and mak- 
ing use of the identity j£ V op t(t')e^ & v„vt(t")dt" df i = 
h (ei St Va Pt (t')dt' _ we find 

X 2 (t) = A 2 (r)e 4(wA(t ~ T) ~^ A 2(*')<it') e -^/TVopt(t')dt' 
+ SlXi(T)e iux( - t - T ' ) e i ^' Ai(t')dt' ei / t T ' A 2 (t')dt' 



x (1- e" 



■f*V opt (t')dt' 



)• 



(15) 



Taking into account the A periodicity Aj(T') = Aj(0), we 
obtain 

A2 ( Q ) = A 2( T ) r -i fj' V mt (t)dt i fj' (A 1 (t)-A 2 (t))dt 

Ai(0) Ai(T) 



— i 



(l_ e -^/r' ' V=Pt (t)dt) > 



(16) 



which, in the limits 



T' 



KptWd* » 1, 



(17a) 



Kpt(t)d*» / [3(A 2 (t))-3(A 1 (i))]^ (17b) 



and for A 2 (T) and Ai(T) of the same order, leads to 

A 2 (0) ^ c, 
Ai(0) ^ ci' 



(18) 



We remark that, denoting the state- vector of the original 

'ai(t) ' 



TDSE |W(t)) 



a 2 (t) 



, the connection to a single 



Floquet vector © leads to X 2 (T)/X 1 (T) = a 2 (T)/a 1 (T), 
i.e. to the ratio of the amplitude at the end of the process. 
If this ratio becomes very large, which corresponds to the 
specific case of an efficient population transfer to state 2, 
the condition (117b[) is not sufficient. It should be replaced 
in general by the condition: 



T' 



T' 



, , V opt {t)dt^> / [S(A 2 (t))-9?(Ax(t))]^ 
n .It jt 

+ log(aa(T))-log(oi(T)). 



(19) 



This is discussed in more detail in Section V.B. 

For an initial condition as a single state of the free 
system, ie. c 2 = 0, c\ = 1, one recovers A 2 (0) <C Ai(0) 
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In this case, we must note that conditions (1TH) are 
less restrictive than condition ©. This is due to the 
fact that conditions (|17j) are obtained constraining a ratio 
of two components, whereas in Q we wished to absorb 
the components, with an error lower than the computer 
accuracy. If the conditions (fl7j) are satisfied, then we can 
force any eigenstate |A) to obey the final condition 



Ai(0) = Al (T) e ^( T '- T )-^' A i« di ) 



(20a) 



A' 



,(0) - A 1 (T)e < ^( T '- T )-J"T T ' Ai(i)d*)£2 (m) 



Cl 



Thus, apart from a global constant Ai(T) which results 
from the diagonalization procedure, an exponentially de- 
creasing term and a global phase, we obtain 



|A(f = 0)) cx m = 0)). 



(21) 



This approximate proportionality is sufficient to impose 
the required initial connection to the Floquet eigenvec- 
tor (HH). This will be illustrated by an example given in 
section IVU 



IV. 



ISOLATING ONE EIGENVALUE IN THE 
FLOQUET SPECTRUM 



The second role of the absorbing potential is to di- 
late the Floquet spectrum and so isolate the "connected" 
eigenvalue hu)\ (i.e. the one associated to the eigenvector 
|A(£)) connected to the initial condition) from the other 
eigenvalue (denoted fruj\> associated to | A' (£))). We con- 
sider for simplicity the initial condition as a single bound 
state |1) of Ho. The absorbing potential takes the form 
set out in Eq. 0. 

We start connecting the solution \ty(q,t)) to the Flo- 
quet vector. This is achieved by solving the stationary 
problem (in the first Brillouin zone): 



for t £ [0, T] : 

d_,(&i(t) n(t) 
at + \ n*(t) A 2 (t) 

for t £ [T, T'] : 

d . f A^t) 



\\(t))=u> x \\(t)), 



(22a) 



dt 







A 2 (t) - iV opt (t) 



|A(i)) = « A |A(t)). 

(22b) 



In the region t £ [T, T'] , we obtain from (|22b|) (see the 
preceding section): 

Ai(i) = Al ( T )^(^(t-r)-^A 1 (t')dO ! (23a) 
A 2 (t) = A 2 (T)e i ^( t - T )-^ A2 ( t ')' it ')e^*^^p t (t')<it'(23b) 

A. Decoupled channels 

The situation is the easiest to follow in the elementary 
case in which the channels are not coupled (fi(t) = 0) and 



with constant diagonal terms Aj. Thus we make the in- 
stant T coincide with t = 0, to study the influence of the 
optical potential alone on the interval [T — 0, T'] without 
any physical coupling terms. In this particular case, with 
T = and t = T' the previous system becomes: 



Ai(T') = A 1 (0)e 4(l ^- Al)T ', 



(24a) 



A 2 (T') = A 2 (0)e 4( ^- A2)T 'e-*-''o r ' fopt(f)«tt\ (24b) 

The same equations can be written for the other eigen- 
state | A'). Floquet eigenvectors must be periodic, i.e. 
Aj(T') = Ai(0). Thus each Floquet eigenvalue must sat- 
isfy simultaneously two conditions: 

1 = e i(cox-Ax)T' ifAl ( )^ (25a) 

1 = e i((^_A 2 )T'+| jf v opt (t)dt) j£ A2(0) ^ o (25b) 

The only solution is to have only one non-zero compo- 
nents for each eigenvector: 

Ai(0)^0 and A 2 (0) = i.e. u\ = A x (26a) 
Ai(0) = and A' 2 (0) ^ 

T' 



2 

i.e. oj\>= A 2 - — y V opt {t)dt 



(26b) 



The terms are not mentioned because we work in 
a given Brillouin zone. In this simpliest case, the exten- 
sion to a N— dimension system is straightforward: all the 
eigenvalues connected to absorbed channels possess an 

imaginary term proportional to J V op t(t)dt. Thus 
we expect to obtain a dispersion of the eigenvalues in 
the complex plane which will leave the other eigenvalues 
distant from the "connected" eigenvalue uj\. 



B. General case 

In the present case of a 2-level coupled system de- 
scribed by Eq. (f22|) . it is possible to go further in the 
analytical description. In the region t £ [0,T], one can 
rewrite (I22al) as 





''Ft 



Ai(t) n(t) 

n*(t) A a (i) 



|A(t))e 



= o, 



(27) 



that is as the same form of the original TDSE of solu- 
tion |\Kt)) = ( a ,l V We connect the two solutions 
V a 2(t) J 

invoking the initial conditions ai(0) = 1, a 2 (0) = 0, and 
Ai(0) = x^Ty^^T'-T)-^' A 1 {t)dt)^ As ( ) _ ( from 
the preceding section): 



ai(t) 

03 (t) 



Ai(T)e 



%{u x {T>-T)-ff Ai(t')rft') _ 



(28) 



The latter equation is just the proof of the Floquet the- 
orem for our specific two-state problem. Considering the 
final physical time t — T, we get 



ai (T) 



(29) 
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that is we connect the imaginary part of the eigenvalue 
lu\ to the final probability amplitude: 

5 (^) = ^7 J 3(Mt))<ft+^log(|oi(T)|). (30) 

To get the counterpart relation for the other ("non- 
connected" ) eigenvalue ujy , we reformulate the complete 
calculation with the adjoint of Hp(t) (using d\ — —dt) ■ 



Hl(t) = Hi + W\i) + V f (i) - ih 

of eigenstates {|Aj j71 (i))} 

Hl.\X j , n {t))=hu>lJ\ j , n {t)), 



d_ 



(31) 



(32) 



where (.)* denotes the complex conjugate. For real en- 
ergies of H and real elements in W(t), this latter equa- 
tion corresponds to the same original problem as before 
but with the use of an exponentially diverging poten- 
tial V'(t). We have then for the components of |A'(f)) 
(denoted as the eigenvector associated to the eigenvalue 
hujy, \X'(i)) is different from \X'(t)) in general): 

X[(t) = Ai(T)e i K'(*- ; D-JjAi(f)«if) > (33a) 

x' 2 (t) = x^TyW'^-ti *z<?w e *£ v '» t M' u 'tf3b) 

which leads in the limits (IT7|) to 
Ai(0)<A^(0), 

A 2 (0) = X' 2 (T)e^l^ T '- T ^lT *2(t)dt) 



o+TlJt V opt (t)dt 



(34a) 
(34b) 



It corresponds to the Schrodinger equation 



a (Am n*(t) 



a' 2 (t) 1 ' 



(35) 



with the initial condition 0^(0) = 0, a' 2 (0) = 1 for which 
we get 

a' 2 {T) = e l tf' A ^ dt e- l ^' T 'e~iJT V ° P *(t) dt . ( 36 ) 

One can connect it to ai(T) as described in appendix lAl 
which induces 



= - l J^[A 1 (t)+A 2 (t)]dt -i ft' A 2 (t)dt 



e +i^ x ,T' e -i-g' V opt (t)dt_ ,^7) 



Identifying (gSJ) and leads to 

i-T' pT 



9f(wv) 



1 
T 7 



3(A 2 (t'))di' + / 9f(Ai (t'JJdf 



7 log(| ai (T)|)-— / V opi (t')df, (38) 



which gives a relation between the imaginary parts of the 
two eigenvalues: 



1 



T' 



V opt {t')dt' - 9f(w A ) 



T' 



T' 



9f(Ai(0 + A 2 (t'))*'- (39) 



This central relation shows that the connected eigenvalue 
will be in general well isolated from the other one for 
a large enough area of the absorbing potential. More 
precisely we have — SJ(wy) 3> — S(wa) when 



1 

KT 



V opt (t')dt' » -29f(w A ) 



T' 



+ - / 3(Ai(t') + A 2 (i'))dt'. (40) 



This feature will be useful in numerical calculations; in 
particular it will improve the rate of convergence of the 
wave operator method Q when applied to the location 
of the thus isolated connected eigenvalue. 

However the separation between the imaginary parts of 
eigenvalues can be not so efficient in practice for specific 
cases of good population transfer. This is analyzed in the 
following section. 



V. NUMERICAL LIMITATIONS AND 
ACCURACY 

In this section we study the numerical limitations of 
the method, restricting the discussion to the situation 
ci(0) = 1,C2(0) = 0. We consider for simplicity the situ- 
ation S(A 2 (t)) = S(Ai(i)) = 0. 



A. General cases 

The accuracy of the method can be roughly estimated 
from the imperfect initial connection with the eigenvector 
| A), that is from the small quantity A2(0). In general, 
when A2(T) and Ai(T) are of the same order, we obtain 
for the error in the final amplitude from (|16|) : 



oi (T) 



(CATM) 



(T) 



4/t V opt (t')dt' 



(41) 



where 



(CATM) 



(T) is the probability amplitude of state 
1 at the end of the physical process obtained from the 
CATM method. This is shown to give a correct esti- 
mation of the accuracy of the method when it is tested 
numerically (see next section). 

We remark that this estimation does not obviously take 
into account the grid size effect. This is studied numeri- 
cally in the next section. 
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B. Case of good population transfer 



A. Some results for selected examples 



The estimation (|4ip is not valid when the population 
transfer at the end of the process is efficient: |ai(T)| — > 0, 
since, in Eq. |Q£), we have then |A 2 (T)/Ai(T)| > f . The 
area of the optical potential should be large enough to 
satisfy the connectivity to a unique Floquet eigenvector: 
A 2 (0)/Ai(0) ~> 0, that is, from {T9| 



The hrst example is a two-state system {|1), |2)} which 
is subjected to a pulsed coupling of frequency little de- 
tuned with the transition frequency. The detuning is 
denoted A and il is the coupling (Rabi frequency). In 
the dressed state picture of the RWA the Hamiltonian is 
(in units such that h = 1) 



T' 



V op t(i)dt»-log(ai(T)). 



(42) 



One limiting case is when there is no separation between 
the imaginary parts of the eigenvalues: 



S(wa') = 9(wa), 



(43) 



leading to 



\£ V opt (t')dt' = -21og(|4 CATM) (T)|)). (44) 

This equation shows that, in this case of equal quasiener- 
gies, the inequality (j42|) is satisfied with only a factor 2. 
More precisely, we have 



A 2 (0) 
Ai(0) 



e ~hSt V opt (t)dt_ 



(45) 



Thus, one can still satisfy A2(0)/Ai(0) — * to get the 
connection to a unique Floquet eigenvector to a good 
accuracy by imposing 



Y h £ v opt (t)dt»i. 



(46) 



This condition (|4l)|) . a bit more restrictive than (|17ap is 
thus sufficient to obtain a quite good relative accuracy 
of the solution in the case of good population transfer, 
even if in that case the imaginary parts of the Floquet 
eigenvalues are close together. 

We can use this limiting case (J43J) to estimate the abso- 
lute accuracy of the method. Assuming $s(u>\i) < Sj(wa), 
we get 



(CATM) 



(T)| > e-™&' v °v^ t ') dt ' i 



(47) 



that is we cannot obtain numerically a population 
|a 1 CA ™^(T)| 2 of state 1 at the end of the physical process 

smaller than ■/? ^optC* ) dt f which gives thus a numer- 
ical limitation of the depopulation of the initial state. 



VI. NUMERICAL INVESTIGATION 

The method is investigated numerically in this section 
through the examples of two- and three- state systems 
driven by a time-dependent field. They can correspond 
for instance to atoms submitted to resonant laser pulses 
in the rotating wave approximation (RWA) [lfl, LUj . 



H = 



ft 

n a 







J7q sin 



. T , 



fi sin 2 (^) A o cos(^ + o 



(48) 

We will consider as initial condition (i) \^(t = 0)) = |1), 
from which we expect a final quasi-inversion of popula- 
tion for larg e enough HqT and AoT (adiabatic passage, 
see [13, and (ii) the more complicated situation 

|*(0)>=ci|l) + c 2 |2). 

The second example is that of a 3 level system 
{|1}, 1 2), 1 3)} driven by two near- resonant laser fields with 
Rabi frequencies Q p and 57 s , tuned to the transitions 
1 o 2 and 2 H 3 respectively. We allow a detuning 
A between the transition frequency 1 — > 2 and the laser 
frequency and assume a two-photon resonance. The ini- 
tial state is |1). Here the RWA Hamiltonian takes the 
form : 



o n p 
ii | n p a n s 

o n., o 



(49) 



We study two situations, on one hand the intuitive case: 
we first turn on the coupling between levels 1 and 2, then 
between levels 2 and 3, 



i2 K 



Tit 



Zi 



fi sin^ — Vt G [0, Ti] (0 elsewhere) 



fl s = Ho sin 
A = A 







) Vt G 


2 L ' 2 



(50) 



With floTi — 20 and AoTi = 0, we expect to observe os- 
cillations without complete population exchange to state 
|3). With A Ti = 20 a partial transfer to |3) occurs with 
less oscillations. 

On the other hand the STIRAP case (Stimulated Ra- 
man Adiabatic Passage) is exactly the inverse of the first 
configuration [l2l |: 



fip = r2p sin 



. 2 /7rf-Ti/2 s 





'1 3 


) Vt G 


.2 U 2 \ 



<t <>::M.r^) V7G[0..T|] 

A = 



(51) 



With fi Ti = 20 and A Ti = 0, STIRAP allows a large 
transfer of the population to |3). 

The total physical time interval T here is 3/2 times 
the period T x of the sine function [0, T] = [0, 3/2Ti]; the 
additional time interval will begin at 3/2T\ for a duration 
of T x . 
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In the subsequent discussion we use the labels (i) and 
(ii) for the 2 level system with initial state |1) and the su- 
perposition of states, respectively. The label (hi) and (iv) 
refer to the 3 level system in the "intuitive" or STIRAP 
situations, respectively. 



B. Calculating with CATM 

From a technical point of view the calculation involves 
the five following steps: 

• Construction of the matrix representation of the 
Floquet Hamiltonian (some details are given in ap- 
pendix |B| 

• Diagonalization of the Floquet matrix 

• Selection of N Floquet eigenstates belonging to the 
first Brillouin zone (for a problem with N levels) 

• Detection of the appropriate "connected" Floquet 
eigenstate, i.e. corresponding to the smallest imag- 
inary part of the eigenvalue in absolute value as a 
criterion 

• Production of the wavefunction via Eq. ([5]) 

In principle only one vector computation is needed. For 
our small-scale examples we can easily use direct com- 
plete diagonalisation. However, for larger systems the 
time-dependent wave operator can be used to find the 
required eigenstate of the corresponding large matrix. 



C. A comparison with direct integration 

We analyse the results obtained with the Floquet 
eigenstate which possesses the smallest value of |3(u;a)|, 
as predicted by the theory. Next we calculate the popu- 
lations 

Pn (t) = \(n\*(t))\ 2 (52) 

and the relative phases 

A»(f)=arg«n|¥(t)» (53) 

for all the previously presented situations. We compare 
the CATM results for the population and phase with 
those of a direct integration using the propagation equa- 
tion 

|tf (t + At)) = e - m ~ lH ( f+ ^r) At \y{t)) (54) 

with At a sufficiently small time-step. For the CATM 
calculation, the size of the Fourier basis set was N = 256 
which is ample for both stable computation and graphical 
representation. 



(a-l) (all) (a-lll) 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(b-l) (b-ll) (b-lll) 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Time t (in units of T) 



FIG. 1. Evolution of populations (a) and phases (b) for the 
2-level system (i) with the initial state |1), and floT = 10 and 
AoT — 10. Exact [i.e. numerical with the direct integration, 
cf. (|54[) ] results (pi and j3\ : short dashes, P2 and /3a : dots) and 
CATM results (pi and /3i: solid line, p2 and fe: long dashes) 
for various amplitudes Vb of the time-dependent absorbing 
potential: (I) VoT = (II) VoT = 10 (III) VoT = 40. 

1. Two-state model 

For the 2 level system (i) the results are shown on Fig. 
[T] In frames (a-I) and (b-I), it is evident that without 
the absorbing potential the use of a single Floquet state 
is not sufficient. On frames (a-II) and (b-II) we can ob- 
serve the effects of the absorbing potential. The initial 
populations approach pi(0) = 1 and p%{Q) = 0, show- 
ing however a small difference of a few percent from the 
reference calculation results. Phases begin to agree with 
those of the reference calculation but the difference re- 
mains important, especially at the beginning. For the 
last case (a-III and b-III), one cannot detect any differ- 
ence between the CATM and the reference results at the 
scale of the figure. 

Fig. [5] shows the same quantities for the initial con- 
dition |*(0)) = ci|l) + c 2 |2), ci = v / 0?75 and c 2 = 
Vfh25- We have used the absorbing potential given by 
Eq. (fTTj) . The previous comments about the efficiency of 
the method remain valid. Fig. [2] illustrates clearly the 
efficiency of the chosen matrix in reproducing the bound- 
ary conditions. 

We now give a more precise analysis of how the exact 
solution is approached. To this end we define a measure 
of the difference between the CATM results and the di- 
rect integration results. For the single component 
calculated by the two methods we define the integrated 
difference of population and of angle: 

*P = ^[ (I(1|*W)catm| 2 - |(1|*W)| 2 ) dt (55a) 

e a = - jf [arg«l|*(t)>cATM)-arg«l|*(t)»]rft(55b) 

These quantities are represented on Fig. [3] as a func- 
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(a-l) 



(II) 
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O 
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(b-ni) 
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0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Time t (in units of T) 

FIG. 2. Same as Fig. [TJ but for the 2-level system (ii), and 
Q.T = 10, AqT = with the initial state 70.75 1) + 70.25 |2), 
for various amplitudes Vb of the time-dependent absorbing 
potential: (I) V T = (II) V T = 10 (III) V T = 40. 
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FIG. 4. Evolution of the populations p n (t) the three- level 
model (iii); exact (numerical) results pi (dots), P2 (long dot- 
dashes), P3 (dot-dashes) and CATM results pi (solid line), 
pi (long dashes), pz (short dashes) without detuning and for 
various amplitude of the absorbing potential, I: VbTi = 0, II: 
5, III: 10, IV: 40. 








. 10 I 1 1 1 1 - 1 1 

10 20 30 40 50 60 

Absorbing potential amplitude V (in units of 1/T) 

FIG. 3. Integrated logarithmic error estimation between di- 
rect integration and CATM with 512 time grid points as a 
function of Vo, calculated with the first component (l|\l/) for 
the 2-level system (i). Errors on population, see Eq. (|55a[) 
(solid line), and on angles, see Eq. (|55b[) (dashed line). We 
remark that these errors follow the anticipated exponential 
law Eq. (|41|l (dotted line) until they reach a plateau due to 
grid effects of CATM 



tion of the absorbing potential amplitude Vq. With 
the logarithmic scale, we observe a quasi-linear law for 
Vq € [10,35] in consistency with Eq. (|4Tj) . The error 
estimates next reach plateaus which are interpreted by 
the grid effects due to the finite basis representation of 
the time in the CATM method. Indeed we can lower the 
level of the plateaus by increasing the number of the grid 
points (not shown). 



(I) (N) 




0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4 



(I") (IV) 




Time t (in units of T) 



FIG. 5. Same as Fig. 0] but with a detuning A Ti = 20. 



2. Three-state model 



For the 3 level system the evolution of the population 
in the three- level model (iii) (as defined in section fVI Aft 
is shown in Fig. |4]and Fig. El without or with detuning 
(A Ti = or A Ti = 20). The selected field amplitude 
was QqTi = 20 and the absorbing potential was gradually 
turned on from VoTi = to VqTi = 40. Here again, if 
the absorption is not sufficient, the results are wrong. 

The result s for the STIRAP model (iv) (as defined in 
section fVI A[) are displayed in Fig [5] The coupling terms 
between levels 2 and 3 are turned on before the coupling 
terms between 1 and 2 and a relatively large population 
inversion is observed. 
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(ii) 
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FIG. 6. Evolution of the populations p n (t) = \{n[$> {t))\ 2 in 
the STIRAP model (iv) ; exact (numerical) results p\ (dots), 
P2 (long dot-dashes), p3 (short dot-dashes) and CATM results 
pi (solid line), P2 (long dashes), P3 (short dashes) without 
detuning and for various amplitude of the absorbing potential, 
I : VqTx = 0, II : 5, III : 10, IV : 40. 
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FIG. 8. 9(cjai) versus Vo for the three- level model (iii). 




10 20 30 40 50 60 70 80 
V (in units of 1/T) 

FIG. 9. Sj(ua 2 ) (solid line) and Q(u)\ 3 ) (dashed line) versus 
Vo for the three-level model (iii). 



10 20 30 40 50 



FIG. 7. (a) (b) Q(uj x ), (c) K(wv) and (d) S(wy), ver- 

sus Vo (all in units of 1/Ti) for the two- level model (i). When 
Vo grows, ujv moves away from uj\ acquiring a imaginary part 
proportional to Vo. 



D. The expansion of the spectrum 

We now analyse the effect of dilatation of the eigen- 
values by the the absorbing potential, that is the feature 
of separating the "connected" eigenvalue with respect to 
the other ones. Fig. [7] shows the Floquet eigenvalues 
{ll>x} and {oj\>} in the first Brillouin zone calculated for 
the two-level model (i) as functions of Vo . Apart for 
small absorbing potential amplitudes where one notices 
an ambiguity concerning the labelling of eigenvalue [l]| , 
3(wa) is a constant value in agreement with Eq. (|30[) . 
and Q(uj\i) shows a linear evolution as predicted by Eq. 

Figures [5] and prefer to the three- level system (iii) and 
show the same features. Concentrating on the imagi- 



nary part of the "connected" Floquet eigenvalue (Fig. |8]) 
we see that, after a region of stabilization, ^(u;;^) is no 
longer affected by the growth of the absorbing potential. 
By contrast both 3(wa 2 ) an d ^{u\ 3 ) acquire imaginary 
parts which are linear with respect to Vq. 

This feature will be useful in practice for large systems, 
in particular if a wave operator method is used to find 
the Floquet eigenstate Q, since that method is efficient 
in finding isolated eigenvalues. 

E. The influence of the number of Fourier basis 
functions 

We give here some details about the stability of the 
results as the number of Fourier basis functions N is 
reduced in the CATM calculation. To increase calcu- 
lational speed and to decrease memory requirement it 
appears necessary to use as small value of N as possible. 
Figure [TU] shows how the final populations obtained in 
the CATM calculation varies as N is increased. These 
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calculations correspond to the STIRAP model (iv). 
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1 0.01 r 

Z> 
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Ll 

0.0001 r 
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N 

FIG. 10. Stability of final populations p\ (solid line), P2 
(dashed line) and P3 (dotted line) functions of the Fourier 
time grid point number N, in logarithmic scale. 

The values p$ ~ 0.82 and p\ ~ 0.18 are stable for 
N > 30 but p-2 ~ 10~ 3 is not obtained accurately until 
N is about 80. We see that finding small probabilities in 
absolute values requires a more precise description of the 
temporal evolution; however about 80 grid points appear 
to be ample for the calculations. The general principle 
is to choose an N which is high enough to follow the 
time variations in the Hamiltonian and to obtain accurate 
values for small probabilities. 



VII. CONCLUSION 

The optimum computational implementation of the 
CATM is still being actively researched but the basic 
principles behind the method are simple to follow. A 
static absorbing potential is often used in treating the 
time development of a wavefunction within the Floquet 
formalism. The novelty of our approach is that the ab- 
sorbing potential is given a time-dependent form such 
that it actively constrains the wavefunction, both by im- 
posing the correct boundary conditions on it and by mod- 
ifying the spectrum so that the specific eigenvalue which 
is appropriate to describe the dynamical process is ren- 
dered relatively isolated from the other eigenvalues. The 
dynamical problem is then rendered into an eigenvalue 
problem in which the isolated eigenvalue is easier to find 
by techniques such as the Bloch wave operator method. 
That it is indeed possible to choose the time-dependent 
potential so as to produce the favourable features de- 
scribed above has been demonstrated for two small-scale 
systems for which accurate comparison results are avail- 
able. For these small test systems the CATM gives ac- 
curate results, although it is clear that the eigenvalue 
problems which arise can involve strongly non-Hermitian 
matrices. 

The CATM has some formal advantages for systems 



with a time-dependent Hamiltonian. A common ap- 
proach for such systems is to use a step-by step propa- 
gation procedure with very small time steps. Many time 
steps are thus required to cover a given time interval and 
this leads to an accumulation of errors as the propaga- 
tion proceeds. By contrast, in the CATM the solution is 
global over the full time interval and so there is no accu- 
mulation of errors; this feature is similar to that shown 
by the (t,t r ) method As expected (and confirmed 
by the present study) the accuracy achievable within the 
CATM is governed by the ability to reproduce the initial 
conditions by suitably adjusting the time-dependent po- 
tential and by the use of a sufficiently dense Fourier time 
grid to describe any fast time variations contained in the 
Hamiltonian. 

Our model calculations have also made clear the role 
of the time-dependent absorbing potential in dilating the 
Floquet spectrum so that the dynamical problem of prop- 
agation within a Hilbert space of a given dimension can 
be converted to that of locating an isolated eigenvalue of 
a non-Hermitian matrix of much larger dimension. The 
difficulty of solving the dynamical problem is thus con- 
verted into the technical problem of devising efficient al- 
gorithms for large non-Hermitian matrices. In Ref. a 
previous version of the CATM was successfully tested on 
a molecular system involving a few hundred states. At 
the moment we believe that the task of isolating and then 
calculating the important dynamically relevant complex 
eigenvalue is probably not possible within the CATM for 
systems which are much larger than those treated in [7|; 
nevertheless, the method may be useful for some systems 
which cause difficulties for the usual propagation meth- 
ods. 
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Appendix A: Properties of dissipative propagators 

We consider a traceless time-dependent dissipative 
Hamiltonian Ht, i.e. having complex diagonal ele- 
ments (with negative imaginary parts) and being self- 
adjoint when it is restricted to its non-diagonal ele- 
ments. It has the corresponding propagator UH T (t,to): 
i-§^UH T (t,t a ) = HTUH T (t,t a ), and its adjoint Hj, the 
propagator U H i (t,t ). They satisfy 

det[Uff T (t,to)] = det[U H i (t,t )] = 1. (Al) 
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^From the definition of the propagators, we obtain 



4[t/L {t,t )U HT (t,t )} = 0, that is 



For the two-state case of general Hamiltonian 

Ai SI 



H 



si* a 2 



(A2) 



(A3) 



with in general complex A x and A 2 and time-dependent 
parameters: Aj = Aj(t), SI = Sl(t), we first decompose it 
as a term proportional to identity and a traceless term: 



Ai + A 2 
H = Z t + H T . 



with 



A 2 -Ai 
2 

si* 



si 

A2-A1 



(A4) 



(A5) 



The propagator for H reads 

CMMo) = e~ l i^o dt '^^ +A ^U HT (t,t ). (A6) 

If Ai and A 2 are real (non dissipative self-adjoint Hamil- 
tonian), the propagator [/# T (i,io) is unitary and is thus 
of the form: 



Uif T (t,to) 



a -b* 
b a* 



(A7) 



In the general case of complex Ai and A 2 , this is not 
true anymore. We write the propagator as 



U HT (t,to) = ( I d J , ad-bc=l. (A8) 
For the adjoint of H: 

gt = AI + AS 1 + g t a (A9) 

the propagator writes 

U H ,(t,t ) =e-*iJ"o t *'(^(0+A3(0)^( t> f ) > (A10) 

where U H t (t, to) connects with Ujj T (t,to) as 



C^t,(t,to) = 



(All) 



These properties are used to obtain a link between 
two wavefunctions resulting from the two orthogonal 

initial slates [ J j and ( ^ j and respectively driven 



by _ff and iff, i.e. UH T (t,t ) ( ?. 



^t(Mo)i ;) = (/ 



« 



and 



Appendix B: Construction of the Floquet 
Hamiltonian 



In this Appendix, we give some details about the struc- 
ture of the Floquet Hamiltonian. For the time dimension, 
we work with a discrete variable representation (DVR) 
{|£i)} ; i = 1 ■ ■ ■ N , associated with a Fourier finite basis 
representation (FBR). We show that the time derivative 
operator can be expressed simply in the DVR basis: Let 
Ij be a column vector of component Sij, i = 1 • • • N, then 



d 



u ) = fft; 



\lu n FFT n (I 3 ) 



(Bl) 



where FFTi represents the i th fast Fourier transform 
component and cjj is the Fourier angular frequency de- 
fined by 



w„ = §5(n-l 



. 1 < ra < f , 
P(n-l-JV) f <n<V. 

Due to the periodicity, (|B2|) is equivalent to 



LJ r . 



U(n-1) l<n<N. 



(B2) 



(B3) 



The matrix representation of 



' dt 



is diagonal in the 



molecular basis, and H(t) and V(t) are approximately 
diagonal in the DVR time basis. Consequently the Flo- 
quet Hamiltonian for the two-level models (i) or (ii) with 
the initial condition c\ = 1, c 2 = is approximately 
represented in the {|1), |2)} ® {\ti)} basis by : 



/ d tll Sl x d tl 
Sl t {d tll + Ai — iVx) 
ft 21 



ft. 



ft 22 
Sl 2 





ft 12 
fi 2 

(ft 22 



*K 2 ) 



\ 



(B4) 



with 



d U 3 
Sli 
A, 

-iVi 



ih 



d_ 

"dt 



Sl(U) Vi,e[0,T] (0 elsewhere) 
A(U) VUe[0,T] (0 elsewhere) 

U - T 



-iV op t(U 



-iVo sin 



T'-T 

\/U e [T, T'\ (0 elsewhere) 



This construction can be directly generalized to treat 
three-level or larger systems. 
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